home *** CD-ROM | disk | FTP | other *** search
/ The Fatted Calf / The Fatted Calf.iso / Applications / Graphics / ContourPlot / Source / splin2.c < prev    next >
C/C++ Source or Header  |  1993-01-26  |  4KB  |  172 lines
  1. /* 2-D bicubic spline from Numerical Recipes in C by Press. W.H et al. */
  2. /* ###### BUG NOTE ######
  3.    NOTE that there is a bug in Numerical Recipes (in both 1-st and 2-nd editions)
  4.    in file splin2.c which is used here with the bug fixed.
  5. */
  6.  
  7. #include <stdio.h>
  8. #include <stdlib.h>
  9. #include "splin2.h"
  10. /* #include "version.h" */
  11.  
  12. #ifndef APPNAME
  13. #define APPNAME "ContourView"
  14. #endif
  15.  
  16. /* float *ytmp,*yytmp;  for debuggin */
  17.  
  18. /* Given x1a, x2a, ya (Data), and y2a (computed by ssplie2),
  19.  * this routine computes interpolated value y at point x1, x2 by
  20.  * bicubic spline.  From Press et al., Num Rec for C.
  21.  * Original in the book is in error.  See the code below.
  22.  */
  23. void splin2(float *x1a, float *x2a, float **ya, float **y2a,
  24.         int m, int n, float x1, float x2, float *y)
  25. {
  26. int j;
  27. float *ytmp,*yytmp;
  28.  
  29.     ytmp=fvector(1, m);    /* was fvector(1, n) which is incorrect */
  30.     yytmp=fvector(1, m);    /* was fvector(1, n) which is incorrect */
  31.  
  32.     for (j=1;j<=m;j++)
  33.         splint(x2a, ya[j], y2a[j], n, x2, &yytmp[j]);
  34.     spline(x1a, yytmp, m, 1.0e30, 1.0e30, ytmp);
  35.     splint(x1a, yytmp, ytmp, m, x1, y);
  36.     free_fvector(yytmp, 1, m);    /* was free_fvector(yytmp, 1, n) which is incorrect */
  37.     free_fvector(ytmp, 1, m);    /* was free_fvector(ytmp, 1, n) which is incorrect */
  38. }
  39.  
  40.  
  41. /* From Press, et al.  Num Rec for C.
  42.  *    Returns second derivatives in array y2a[1..m][1..n]
  43.  */
  44. void splie2(float *x1a, float *x2a, float **ya, int m, int n, float **y2a)
  45. {
  46. int j;
  47.     for (j=1;j<=m;j++)
  48.         spline(x2a,ya[j],n,1.0e30,1.0e30,y2a[j]);
  49. }
  50.  
  51.  
  52.  
  53.  
  54. void splint(float *xa, float *ya, float *y2a, int n, float x, float *y)
  55. {
  56. int klo,khi,k;
  57. float h, b, a;
  58.  
  59.     klo=1;
  60.     khi=n;
  61.     while (khi-klo > 1) {
  62.         k=(khi+klo) >> 1;
  63.         if (xa[k] > x) khi=k;
  64.         else klo=k;
  65.     }
  66.     h=xa[khi]-xa[klo];
  67.     if (h == 0.0) {
  68.         fprintf(stderr,"%s: Bad XA input to routine splint()\n", APPNAME);
  69.         exit(1);
  70.     }
  71.     a=(xa[khi]-x)/h;
  72.     b=(x-xa[klo])/h;
  73.     *y = a*ya[klo] + b*ya[khi] + ((a*a*a-a)*y2a[klo] + (b*b*b-b)*y2a[khi]) * (h*h)/6.0;
  74. }
  75.  
  76.  
  77.  
  78. void spline(float *x, float *y, int n, float yp1, float ypn, float *y2)
  79. {
  80. int i,k;
  81. float p,qn,sig,un,*u;
  82.  
  83.     u=fvector(1,n-1);
  84.     if (yp1 > 0.99e30)
  85.         y2[1]=u[1]=0.0;
  86.     else {
  87.         y2[1] = -0.5;
  88.         u[1]=(3.0/(x[2]-x[1]))*((y[2]-y[1])/(x[2]-x[1])-yp1);
  89.     }
  90.     for (i=2;i<=n-1;i++) {
  91.         sig=(x[i]-x[i-1])/(x[i+1]-x[i-1]);
  92.         p=sig*y2[i-1]+2.0;
  93.         y2[i]=(sig-1.0)/p;
  94.         u[i]=(y[i+1]-y[i])/(x[i+1]-x[i]) - (y[i]-y[i-1])/(x[i]-x[i-1]);
  95.         u[i]=(6.0*u[i]/(x[i+1]-x[i-1])-sig*u[i-1])/p;
  96.     }
  97.     if (ypn > 0.99e30)
  98.         qn=un=0.0;
  99.     else {
  100.         qn=0.5;
  101.         un=(3.0/(x[n]-x[n-1]))*(ypn-(y[n]-y[n-1])/(x[n]-x[n-1]));
  102.     }
  103.     y2[n]=(un-qn*u[n-1])/(qn*y2[n-1]+1.0);
  104.     for (k=n-1;k>=1;k--)
  105.         y2[k]=y2[k]*y2[k+1]+u[k];
  106.     free_fvector(u,1,n-1);
  107. }
  108.  
  109.  
  110. /* memory allocation routines */
  111.  
  112. /* ###### 1-d float */
  113.  
  114. float *fvector(int nl, int nh)
  115. {
  116.     float *v;
  117.     v=(float *)malloc((unsigned)(nh-nl+1)*sizeof(float));
  118.     if(!v) fprintf(stderr, "allocation failure in fvector()\n");
  119.     return(v-nl);
  120. }
  121.  
  122.  
  123. void free_fvector(float *v, int nl, int nh)
  124. {
  125.     free((void *)(v+nl));
  126. }
  127.  
  128.  
  129. /* ##### 2-d float */
  130.  
  131. float **fmatrix(int nrl, int nrh, int ncl, int nch)
  132. {
  133.    int i=2,j;
  134.    int error=0;
  135.    float **m;
  136.     m=(float **) malloc((unsigned) (nrh-nrl+1)*sizeof(float *));
  137.     if (!m) fprintf(stderr, "allocation failure 1 in fmatrix()\n");
  138.     else{
  139.        m -= nrl;
  140.        i=nrl;
  141.        do{
  142.         m[i]=(float *) malloc((unsigned) (nch-ncl+1)*sizeof(float));
  143.         if (!m[i]) {
  144.             fprintf(stderr, "allocation failure 2 in fmatrix()\n");
  145.             error=1;
  146.         }else
  147.             m[i] -= ncl;
  148.         i++;
  149.            }while( (i<=nrh) && (!error));
  150.     }
  151.     if(error){
  152.         fprintf(stderr, "Freeing allocated submatrices...\n");
  153.         for(j=i-2;j>=nrl;j--)           /* free those allocated */
  154.             free((void *)(m[j]+ncl));
  155.         free((void *)(m+nrl));
  156.         m=NULL;
  157.      }
  158.     return m;
  159. }
  160.  
  161.  
  162. void free_fmatrix(float **m, int nrl, int nrh, int ncl, int nch)
  163. {
  164.     int i;
  165.     for(i=nrh;i>=nrl;i--)free((void *)(m[i]+ncl));
  166.     free((void *) (m+nrl));
  167.  
  168. }
  169.  
  170.  
  171.  
  172.